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INTRODUCTION 


The classical Kepler equations for the elliptic and hyperbolic 
cases of two-body motion are written 

M - E - e sin E 
M = e sinh H - H 


respectively, where 

M = (t ' t) • 

Here a is the semi-major axis; t is the time of pericenter passage; 
p, = + 1112 ), where G is the universal gravitational constant and 

and ms are the masses of the two bodies; and e is the eccentricity 
of the relative orbital motion. The quantities E and H are defined by 
the relations 

tan E/2 = J tan 0/2 

tanh H /2 = J Y~ tan 9/2 

where 0 is the true anomaly. Note that we have used the sign conven- 
tion for the semimajor axis, i.e., a > 0 for an ellipse and a < 0 for 
a hyperbola. In both cases a is obtained from the equation 



where r and v are the magnitudes of the position and velocity vec- 


tors respectively. 


It is easily shown that Kepler's equation possesses a unique 

solution. When the eccentricity is not close to one, there are effi- 

) 

cient iterative methods for obtaining this solution. However, when a 
becomes large and forces the eccentricity near unity, both the elliptic 
and hyperbolic equations suffer a critical loss of accuracy. These are 
designated "nearly- parabolic" orbits and their solution requires special 
treatment . 

The parabolic motion is described by the cubic equation 


2 ./ ^5 (t - t) =* tan 9/2 + ^ tan 3 9/2 


(5) 


where p is the semilatus rectum and is such that 


p = a ( 1-e 2 ) 


for elliptic and hyperbolic motion and 


P = 2q 


for parabolic motion, where q is the distance between the focus and 
vertex of the parabola. It is easily shown that one and only one real 
root exists for this equation. 


THE GAUSS METHOD 

The classical method for solving Kepler's equation in the nearly 

parabolic region is an iteration approach due to Gauss, and a detailed 

[5] 

description of the technique has been given by Herget . The method 
requires the use of auxiliary functional quantities which are usually 


2 


[21 

obtained from special tables . However, Benima et al L J have derived 
series expansions for these quantities which makes the method suitable 
for high-speed computer solution. Following the notation of Benima, the 
equations for Gauss' solution are 

a -itiie b e - /msr ( 6 ) 

V 10 “l + 9e V l+9e 


A = b tan 2 w/2 


( 7 ) 


B a 


(t - t) 


tan w/2 



w/2 


( 8 ) 


tan 9/2 = c C tan w/2 


( 9 ) 


where 


A = 15(E - sin E) 
9E + sin E 


_ 20 f+K 
9E + sin E 


C = — tan E/2 

JTK 


( 10 ) 


for elliptic orbits and 


A _ 19(H - sinh H) 
" 9H + sin H 


B = 20 /~A 

9H + sinh H 


( 11 ) 


C = -i— tanh H/2 


for hyperbolic orbits . 

Since B and C are functions of A, they may be tabulated or, as 
developed by Benima et al, expanded as 
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00 


00 


(12) 



J=0 J=0 


This form is more efficient for computers, because now only the coef- 
ficients need be stored. The first eight coefficients for these expan- 
sions have been presented by the above authors . 

The Gauss procedure to determine 0 from (t - t) is to solve equa- 
tion (8) for tan w/2 by successive approximation, beginning with B = 1. 
The value of tan w/2 obtained by solving the cubic equation with B = 1 
permits the computation of A by equation ( 7 ) • This in turn yields a 
new value for B by the series expansion (12) and permits a more ac- 
curate solution for tan w/2. This process is repeated until A reaches 
a desired accuracy, and then tan 9/2 is computed from equation ( 9 )* 
having used the expansion (12) for the calculation of C. Rapid con- 
vergence of the method results from the condition that = 0 in equa- 
tions (i2) . Numerical efficiency is increased if, after the first 
step through the algorithm, the cubic equation is solved by a single 
iteration using tne solution of the cubic from the previous step as an 
initial guess. It should be noted tnat double precision accuracy re- 
quirements for other tnan short time intervals would force the calcu- 
lation of more series coefficients than the eight given and hence the 
expense of considerable labor. 

The calculation of (t - t) when given 9, the reverse of the above 
problem, is not rapidly convergent using the Gauss method. For this 


k 


calculation Beniraa et al introduce the series 


C " 2 = I ^ (13) 

J=0 

which is somewhat more rapidly convergent than the series for C. The 
procedure is then to set C =1, obtain tan w/2 from equation ( 9 ) , com- 
pute a value for A from equation (7)> and then calculate a new C 2 
from the series (13) . This scheme is repeated until A reaches its 
final value, then B is computed from equation (12) and (t - t) from 
equation ( 8) . 

The position and velocity may be written 


r = q D (l - tan 2 9/2) 1 F + 2q D tan 9/2 1^ 


(14) 


= /^iTST (l + tan^ 9 ' /g ) L‘ 2 tan 9/2 V U+e • ta “ 2 9/2)I n] 


where is a unit vector in the direction of pericenter and 1 is a 
unit vector in the direction of l x 1^., where L is the angular momentum 
of the orbit . Here 


for an ellipse and 


D = | (1 + cos E) 

D = 7; (l + cosh H) 


for a hyperbola. As before D may be tabulated as a function of A, 
expanded as 



THE COEFFICIENTS 0 . , y , , a , , AND 6 
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♦From Benima, et al 



or simply calculated from 


D ■ l A ' 


D (1 + AC 2 ) 


The coefficients 0 , y . , < 7 . and 6 . are presented in Table II of 
J J J J 

reference [2] and are reproduced in our Table I for completeness. 


THE UNIVERSAL METHOD 

In place of tne separate equations given above to describe ellip- 
tic, parabolic and hyperbolic motion, a universally valid Kepler equa- 
tion may be written^- ^ 

(t-t 0 ) * a b 3 s(<*3 2 ) + b p 2 c(ae 2 ) + r 0 a 


r 0 * v 0 


A = 1 - a tq, 


1 2 v o 

a ~ r 0 p, 


where r 0 and v 0 are the position and velocity vectors at time t = t 0 
and the C and S f notions are defined by the series 


S(x) = I 


(2i + 3) : 


C (x) = y — 
jTq (21 + 2) ! 


The symbols A, B, C and 3 defined above are independent of these of the 
previous section. Given the time t and the position and velocity at 
t 0 , 3 is found from equation ( IT) • The position and velocity at time 
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t are 


n 


? = [i - c (oB 2 )] r 0 ♦ [(t - t 0 ) 
v = £arB 3 S(aB 2 ) “ Bj r 0 + “ 



CCcrB 2 ) 1 v 0 
r _j 


r = A B 2 C(<yB 2 ) + 


B [_8 


■ a B' 


S(crB 2 ) 


r 0 


(19) 


The universal Kepler equation, which is independent of the peri- 
center location has several important features. For a = 0, the para- 
bolic case, equation (17) reduces as it should to the cubic equation 


f(B) = JZ (t - t 0 ) (20 

where 

f ( B) * ^ A p 3 + j B e 2 < r 0 B. (21) 

For the parabolic case 0 may then be identified as /p (tan 0/2 - tan 0o/2) . 
Moreover if we define 


y( B) = A B 3 S( qtB 2 ) + B B 2 C(a0 2 ) + r 0 0 (22) 

we find that 

^ = r 0 ♦ A 0 2 C(aB 2 ) + B [b - a- B 3 S(c*0 2 )] - r (23) 

when is evaluated a-c the root of y( B) = /\i (t-to) • We also find 
ap 

that y(B) has inflection points given by 
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P. „ = — tan ( - ZLl ) 
inf sz ' a / 


when <* > o and 


B inf “ tanh '’ (- A ) 


when q- < 0. Thus there is only one inflection point for hyperbolic 
motion and an infinite number given by 


B inf ' 6 inf + ^ - 
° Jot 


n = G, ±1 , ±2, 


for elliptic motion, where B. r denotes the principal value of the 

ini o 

arctan in equation (24) . The value of y at an inflection point is 

&inf + B 

given by y( 8 inf ) = for all a . 

Various tecnniques exist in the literature for the solution of 
equation (lj) for B when (t-t 0 ) is given. In general an iterative 
method will be found most efficient, but in certain instances other 
methods of solution may be advantageous 1 - J . 

For the case under consideration, that of nearly parabolic motion, 
a is a small quantity. This fact has been utilized in expanding the 

r 3i 

solution of equation (17) as a power series in cr • However, under 
this situation the universal Kepler equation may be solved very effi- 
ciently for 3, given (t-t 0 ) by calculating an initial value for 3 from 
the cubic equation (20) and then applying a Newt on-Raphson iteration 
technique . Optimal computing forms for the C and S functions have been 
obtained^ which greatly increase the numberial efficiency when 
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solving the universal formulation. This method of solution then has 
advantage over the Gauss method in that it is not tied to periapsis 
location, and moreover it is a special solution for nearly parabolic 
motion only in the selection of an appropriate starting value for the 
iteration procedure . In direct numerical comparison with the Gauss 
procedure of Benima et. al, the methods displayed equality for short time 
intervals ( approxinBtely three iterations for double precision accuracy), 
but otherwise the universal technique displayed marked superiority. 
ric”res 1 and 2 are graphs of y( p) and f( p) for elliptic and hyperbolic 
orbits respectively showing the approximating features of the cubics 
f(p) expanded about inflection points. 

In contrast to the Gauss method, the reverse calculation using 
the universal Kepler equation is as efficient as the procedure outlined 
above. If P is a given quantity, (t-to) may be simply evaluated from 
equation (17) with no numerical difficulty. The calculation with 
( 0 - 0 o ) as a given quantity is a little more involved. The dependence 
of the eccentric anomaly and its hyperbolic counterpart on the true 
anomaly may be written 

E-Eo 

tan — — = Jot d 
H-Ho 

tanh — — = ,/ -oi d (27) 

e-e 0 

r 0 sin -r — 

0-0o e-e 0 

/p" cos — — - B sin —5 — 

_ rg 

cot [Vro (2-0*0) - J - B 
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y ( P) = Ap 3 S(«|3 2 ) + 

B0 2 C(«P 2 ) +r 0 P 

f(P)={ A (P -g ;nf ) 3 + 

f(P -P !nf ) 2 + 

r 0 (P ‘ 0 inf> +yOinf) 



II 


where we express 


p = a( 1-e 2 ) = 2r 0 - 0 /r o 2 - B 2 

to retain significant digits which would otherwise be lost when e is 
close to unity. 

Define 

X = , x = -2^2 (28) 

2/oi 2f-at 


for elliptic and hyperbolic motion respectively. There then results 
the universal equation 


x - a x 3 S(cyx 2 ) =* d £l - qx 2 C(qoc 2 )J 


(29) 


where x is related to the universal variable p of equation (17) by 



(30) 


The solution of equation (29) may be easily obtained by a Newton-Rapnson 
method, yielding 


x 


i+1 


x. 

1 


|~x i - a x i 3 s( a x i 2 )j - d£l-Q'X i 2 c(ax 1 2 )j 
[l-aXj 2 cCerx^Jj + adj^Xj-aXj 3 SCax^) j 


(31) 


The authors have found that full double precision accuracy is usually 
attained in only two iterations with the use of 
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( 32 ) 


x 0 = — tan 1 (/a d) 

/a 

xq = — tanh 1 d) 

/=5 

as initial values for the elliptic and hyperbolic cases respectively. 

Then by equations (30), ( 17) and reduction formulas for the C and S 
functions we may calculate the time from 

—■ (t-t 0 ) = 4 Ax 3 S(tarx 2 ) + 2 Bx 2 C(4<?x 2 ) + r 0 x 
or (33) 

^ (t-t 0 ) = Ax 3 S(q'X 2 ) + r 0 x + ^1-cyx 2 S(<*x 2 )J 

'J\x 3 C(cyx 2 ) + Bx 2 ^1-qoc 2 S(ax 2 )^j . 


It should be noted that another technique for the solution of the 

reverse calculation which is more efficient than the Gauss method lies 

in the use of the unified form of Lambert's theorem, due to Iancaster 

[7] 

and Blanchard , where a series expansion is given for the normalized 
time of flight in the near parabolic region. 
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